Breakdown of the Thomas-Fermi approximation for polarized Fermi gases 
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We use Bogoliubov de-Gennes theory to show that the commonly used Thomas-Fermi approxi- 
mation (TFA) can fail in describing polarized unitary gases in anisotropic harmonic traps. We find 
a magnetized superfluid region inside the trap, with order parameter oscillations, even though there 
is no such stable bulk phase. This leads to magnetization profiles that deviate from contours of 
constant potential energy. We determine how this violation scales with trap anisotropy and number 
of particles, and show that we are able to account for important differences between the MIT and 
Rice experiments. 



The study of strongly interacting Fermi systems using 
ultracold atomic gases has attracted interest across the 
physics community. The phenomena observed in these di- 
lute gases are expected to shed light on systems as diverse 
as high Tc superconductors, quark-gluon plasmas and 
quantum chromodynamics. The great virtue of atomic 
gases is the tunability of interactions using a Feshbach 
resonance so that the entire BCS to BEG crossover can 
be studied, with the most strongly interacting, unitary 
regime in-between these two limits. A series of beauti- 
ful experiments have probed condensation of fermionic 
pairs [H, pairing gaps [l], quantized vortices and the 
thermodynamics [J] of the crossover. A particularly ex- 
citing new direction is the study of partially polarized 
gases [HI, [gI with an imbalance in the number of up and 
down 'spin' fermions. 

An important aspect of cold atom experiments is the 
presence of a harmonic trap. A "local density" or 
Thomas- Fermi approximation (TFA) 0] has usually been 
adequate to take this into account. The TFA asserts that 
the properties at point r in a 'slowly varying' potential 
are the same as those of the uniform gas at a chemical po- 
tential /icr(r) = fia — with a =T, j. This leads to the 
simple result that the spatial dependence of any observ- 
able must follow contours of constant trapping potential, 
which is directly testable for the densities ncr(r). 

The two experiments on polarized, unitary Fermi gases 
find rather different results with respect to the TFA. The 
MIT group [I, [11, with a large number of atoms N = 10^ 
and a small trap anisotropy 1/a = ujr/^z — 5, finds 
that the densities follow contours of V{r). On the other 
hand, the Rice group 0, with smaller N = 10^ and 
larger anisotropy 1/a :^ 50 observes gross violations of 
the equipotential contour condition for the densities [loj . 

Motivated by this, we have investigated the validity of 
the TFA using T = Bogoliubov-deGennes (BdG) cal- 
culations fn\ and scaling arguments for the unitary gas 
in anisotropic, three dimensional traps with polarization 
up to 40%. Our main results are: 

(1) The TFA is always violated in a trap in so far as 
the spatial variation of the order parameter is concerned. 
Between the unpolarized superfluid at the center and the 
fully polarized normal gas at the edges of the trap, there 



is an intermediate region which is a magnetized super- 
fluid with an FFLO-like oscillation [12[ of the order pa- 
rameter. We emphasize that there is no corresponding 
stable phase for the uniform system. 

(2) We find that the size of the magnetized superfluid 
region depends on both polarization and anisotropy and 
can be much larger than kj^ for large anisotropy. 

(3) The violation of the equipotential contour criterion 
for the magnetization m(r) = n^{r) — n|(r) increases 
with increasing anisotropy but decreases with in- 
creasing total number of particles N = -\- N^. 

(4) We derive a simple condition for the consistency of 
the TFA: Sj^/cOr = {Naf^^f{P) > 1, where / is a func- 
tion of the polarization P = {N^ — N^) /N. We use this 
{Na)^^^ scaling of dji/ujr to get a better understanding 
of the N and a dependences of our BdG results. 

(5) We are thus able to account for the differences be- 
tween the MIT and Rice experiments with respect to the 
question of when the magnetization should or should not 
follow contours of constant potential at T = 0. 

(6) A general implication of our results is that the TFA 
based on the bulk phase diagram must be used with cau- 
tion to describe a system in a trap when there are several 
competing phases. 

Bogoliubov-deGennes equations: Our approach 
to the problem of strongly interacting, polarized Fermi 
gases in anisotropic traps is to solve the Bogoliubov- 
deGennes (BdG) equations 0. This is the simplest 
approach which goes beyond the TFA and is a gener- 
alization of the BGS-Leggett mean field theory for a spa- 
tially inhomogeneous gas. This method has been ap- 
plied to the study of vortices in the strongly interact- 
ing regime [M], as well as the study of polarized gases 
in isotropic traps For a single-channel description 

valid for the experimentally-relevant wide resonance, the 
Hamiltonian density for the polarized gas is 

H{v) = ^ft(r)[jf„(r)-^^]xl/^(r) 

cr 

-5*|(r)*|(r)*i(r)*T(r), (1) 

where Hq = — V^/2m + V{r), m is the fermion mass 
and we set = 1. The trapping potential is V{y) = 
\muoQ{r'^ ^ o?z^\ where we use cylindrical coordinates 
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FIG. 1: (A) Majority (black) and minority (gray) density 
profiles and (B) order parameter A (black) and magnetization 
(grey) profiles along the axis of the trap. The solid lines 
are BdG results and the dashed lines are TFA results. The 
calculations are for a trap with a = 1/4 containing N = 865 
particles and a polarization of 30%. 



r = {r^O^z). We define the average chemical potential 
jii = (/i| + /i^ ) /2 and the difference as 2/i = = /i^ — /i| . 
The mean field state is found through the solution of the 
BdG equations [lH 



Ho{r)-fi A(r) \ f u,{r) 
A*(r) -i7o(r)+/i ; U.(r) 



Vi{r) 



(2) 

together with the gap equation, polarization and total 
density at zero temperature given by 



A(r) = 9 ^ ^.,(r)<(r), (3) 

Ei>h 

m(r) = Yl {\ui{r)\' + \v,{r)\^), (4) 

0<Ei<h 

n(r) = m(r)+ ^ 2|vi(r)p. (5) 

Ei>h 

These equations are solved self-consistently for A(r),/i 
and h using the constraints that the total number of 
particles is N = J d^rn{r) and the polarization P = 
N-^ J d^rm{r). 

The solution of these equations is simplified if we ex- 
pand the wavefunctions in terms of the eigenfunctions of 
the diagonal piece Ho{t) — /i. Measuring lengths in units 
of the radial harmonic oscillator length ao = 1 / y/mcoo 
and energies in units of cjq, these functions are (j)np£ = 
fp£{r) ex.p{i£6)gn{z)/\^, where the radial and axial 
functions are related to associated Laguerre and Hermite 
polynomials, fp£{r) = ^/p\/Jp~T~^.e~^ /^rl^lL^(r^) and 

gn{z) = v'v^/(2^v^n!)e-"^'/2^n(v^^), respectively. 
The corresponding eigenvalue is enpi = {2p + £ + !) + 



in i due to axial symmetry. For a given 
diagonalize 



aw -t(^) 



we need to 



(6) 



where T^^\ , 

nn'pp' 



^nl'pp'= / I dzfp£{r)fp^£{r)gn{z)gn^{z)A{r,z) 

Jo J-oo 

Since A(r, z) = A(r, — z), the only non-zero matrix ele- 



^nn'^pp' S^nd 



ments of A 



correspond to even n -\- n' . 



a{n + 1/2) — /i. The BdG Hamiltonian is block-diagonal 



nn'pp' 

The bare coupling g in eqs. (|ll3p is related to the two- 
body s-wave scattering length through m/(47ras) = 

1/g + 771^/2^/(^2^2)^ ^^^^^ 

is an ultraviolet en- 
ergy cutoff. The number of particles determines = 

{6Na)^^^uJo and we have used a cutoff of Ec = 4ej for 
the calculations in this paper. 

Densities and Order Parameter: We now dis- 
cuss the self-consistent solution of the BdG equations 
as a function of total N = -\- N^^ polarization 
P = {N^ ~ ^i) and trap anisotropy 1/a at unitarity 
{as = oo). We have extensively studied the problem for 
N up to 2500 particles, < P < 0.4 and a = 1, 1/2, 1/4. 
In Fig. 1(A) we plot the majority (n^) and minority 
(77|) densities along the z axis for a representative data 
set (a = 1/4, TV = 865 and polarization P = 30%). 
In Fig. 1(B) we plot the corresponding magnetization 
777 (r) = n| — 77| together with the local order parame- 
ter A(r). In both panels the solid lines are BdG results, 
while dashed lines are TFA predictions (using the bulk 
phase diagram [16.] as input). 

Both the BdG and TFA results show an unpolarized 
superfluid at the center of the trap and a fully polarized 
normal gas at the edge. There is a marked decrease in 
the BdG central density relative to TFA, with a redistri- 
bution of minority atoms to an intermediate region. The 
main difference between BdG and TFA is precisely in this 
intermediate region. Within the TFA there is a discon- 
tinuous jump in the order parameter which is smoothed 
out in the BdG solution since this lowers the gradient 
energy. What is perhaps unexpected is that the decay- 
ing BdG order parameter exhibits oscillations, similar to 
those expected in a putative FFLO phase [12] with a pe- 
riod roughly consistent with 27r/ {kp^ — kpi) (where kpa 
are the local Fermi wavevectors). Irrespective of the im- 
portance (or otherwise) of the limited number of rather 
small amplitude oscillations, it is unambiguous that this 
intermediate region is a magnetized superfluid: it has 
both a non-vanishing superfluid order parameter and a 
non-zero magnetization. 

We note that the magnetized region that we see is com- 
pletely different from a partially polarized normal region. 
The latter, though barely visible, is present in the TFA 
results of Fig. 1(A) and derives from the partially polar- 
ized normal phase which exists for a small range of jii/ h in 




FIG. 2: Top row: False color plots of the 3D magnetization as 
The trapping potential has a — 1, 1/2, 1/4 for the left, center 
column integrated magnetization densities. 

the bulk mean field phase diagram [iG^ . Such a partially 
polarized normal region, which is also seen prominently 
in phenomenological implementations of the TFA [13] , is 
never observed in our BdG calculations with P < 40%. 

The existence of a magnetized superfluid region im- 
plies the breakdown of TFA. There is no such phase for a 
uniform gas in the thermodynamic limit at unitarity and 
this region is stabilized only by the presence of a trap. 
Earlier BdG studies in isotropic traps [11] already found 
such an FFLO-like region but, as we discuss next, it has 
a much larger spatial extent in the anisotropic case. 

It is important to ask if the intermediate region is suf- 
ficiently narrow that it can be described as an interface 
with a surface energy [IS]. While this may be a rea- 
sonable semi-phenomenological description, we believe it 
is not microscopically correct. From our BdG calcula- 
tions we find that the size of the intermediate region is 
proportional to k^^ but with a large proportionality con- 
stant which is P and a dependent, and increases rapidly 
with trap anisotropy along the axial direction. For 
instance, in Fig. 2 the 2;-extent of the magnetized super- 
fluid is 6 times, 9 times and 16 times the local kp^ for 
a = 1, 1/2 and 1/4 respectively [loj. 

We emphasize that, despite an apparently widespread 
feeling to the contrary, the size of the intermediate re- 
gion is also quite large in the experiments. This size can 
be determined from the separation between the point at 
which and flrst deviate from each other to the point 
where vanishes. Even the MIT data, which shows lit- 
tle evidence for breakdown of TFA, has an intermediate 
region of 10/im along the radial direction, a signifl- 
cant fraction of the unpolarized core radius; see Fig. 5(b) 
of ref. [1]. In the Rice data the intermediate region is 
:^ 100/im along the axial direction. This is two-orders of 
magnitude larger than kj^ and a signiflcant fraction of 
the superfluid core size; see Fig. 3(b) of ref. [9]. 

Equipotential contours: At present there are no di- 
rect experimental probes of the magnetized superfluid re- 
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a function of {z,r) for a system with N = 865 and P = 30%. 
and right panels respectively. Bottom row: the corresponding 



gion. Thus we must look for signatures of the breakdown 
of TFA in the testable question of whether the densities 
follow contours of constant potential. 

We flrst show how the equipotential contour condi- 
tion is progressively violated as a function of increas- 
ing trap anisotropy 1/a. In the top panel of Fig. 2 
we plot the magnetization m(r, z) for N = 865 parti- 
cles with P = 30% polarization and a = 1,1/2,1/4. 
In the lower panel we show the corresponding plots of 
the column-integrated magnetization, which is simpler 
to measure in experiments, and is given by nicoiiy^z) = 
J^oQ dx -\- y'^^ z). For the spherical case {a = 1) 

the equipotential contour condition must be satisfled by 
symmetry, despite the violation of TFA in the order pa- 
rameter. As the anisotropy increases we see that the 
magnetization gets more concentrated along the wings. 
Moreover, the boundary between the magnetized and 
unmagnetized regions near the z = plane becomes 
straighter, yielding a magnetization "hole" that becomes 
more rectangular. This is very similar to the observed 
proflles in the Rice experiments [9]. 

To understand why the MIT results look so different, 
we must study the dependence on the total number of 
particles, for a flxed trap anisotropy and polarization. In 
Fig. 3 we plot the results for (from top to bottom) = 
865, 1538, and 2307 particles in a trap with a = 1/2 and 
P = 30%. For the smallest N the magnetization is local- 
ized along the axis, and is seen to spread out toward the 
radial direction with increasing N. The largest N results 
show a rather elliptical magnetization density indicating 
that the magnetization begins to follow the equipotential 
contours as N increases. 

Scaling with N and a: We next present a simple ar- 
gument for the consistency of the TFA which will allow us 
to see better how our results scale with N and a. In a po- 
larized Fermi gas the conditions for the consistency of the 
TFA are jUa/^o ^ 1 and dfi/uoQ ^ 1. The last inequality 
is the one that is most easily violated when P is not close 
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FIG. 3: Magnetization profiles as a function of total N (from 
top to bottom: 865, 1538 and 2307) at a fixed trap anisotropy 
(a= 1/2) and P = 30%. 

to unity. We start with the TFA densities obtained from 
the mean- field phase diagram at unitarity [16]. By spa- 
tially integrating these density profiles we find and 
Ni in terms of Sja and the average /i. Inverting these 
relations we can show that = {Naf^^f{P) 

Here /(P) is a monotonically increasing function of P 
which goes like P^/^ for P <C 1 and is of order unity for 
the P values of interest. 

We have checked that our BdG results for Sjii/coo are 
consistent with the (Na)^^^ scaling, even though the val- 
ues are smaller than the TFA estimates. If the condition 
(5/i/cc;o ^ 1 is violated the TFA will breakdown. How this 
breakdown manifests itself in the size of the magnetized 
superfiuid and in the violation of the equipotential con- 
tour criterion depends strongly on the anisotropy 1/a, 
as seen above. For a given (Na)^^^ these violations are 
more pronounced for larger 1/a. 

The Rice experiments with (Na)^^^ 10 and a large 
anisotropy show a significant violation of the equipoten- 
tial contour criterion, consistent with our findings. On 
the other hand, we also understand why the MIT exper- 
iments, with (Na)^^^ ^ 100 and a small anisotropy, find 
that this criterion is obeyed. 

In summary, we have shown using BdG equations for 
the unitary, polarized Fermi gas that the often used TFA 
breaks down with the appearance of an intermediate 
magnetized superfiuid region. The deviations from TFA 
are more pronounced with increasing trap asymmetry, 
both in terms of the size of the intermediate region and 
in the violations of the equipotential contour criterion. 
These violations become less important for large parti- 
cle numbers and more spherical traps, in a way that is 



consistent with current experimental results. Important 
questions for future work are inclusion of finite tempera- 
ture effects, ways to probe the intermediate magnetized 
superfiuid region, and understanding the nature of the 
ground state at large polarization. 
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for help with the preparation of the figures. 
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